%%
% 本程序求解静电场问题
%%
close,clc,clear;
% 材料参数设置
epsilon.air = 1;  %---介电常数
epsilon.cablein = 1;
epsilon.cableout = 1;
%% 
% 激励大小设置(V)
voltage_cablein = 1;
voltage_cableout = 0;
%% 
% 读取网格信息
data = load("mesh\T2.mat");  %---导入GmshInf网格文件
coord = data.coord;                     %---坐标信息
coord(:, 2:end) = coord(:, 2:end)*1e-3;  %---坐标单位变换为国际单位制：m
element_v = data.element_V;
element_s = data.element_S;             %---获取面单元信息
element_l = data.element_L;             %---获取线单元信息

% 常数设定
PosOfcoordx = 2;    %--结点x坐标列
PosOfcoordy = 3;    %--结点y坐标列
PosOfcoordz = 4;    %--结点z坐标列
PosOfcoord = 1;     %---结点编号列
PosOfelement = 2;   %---单元编号列
PosOfLentity = 1;   %---线单元实体编号列
PosOfSentity = 1;   %---面单元实体编号列
PosOfVentity = 1;   %---体单元实体编号列
PosVEleNum1 = 3;     %---单元中结点局部编号1对应的结点全局编号列
PosVEleNum2 = 4;     %---单元中结点局部编号2对应的结点全局编号列
PosVEleNum3 = 5;     %---单元中结点局部编号3对应的结点全局编号列
PosVEleNum4 = 6;     %---单元中结点局部编号4对应的结点全局编号列
PosSEleNum1 = 3;     %---单元中结点局部编号1对应的结点全局编号列
PosSEleNum2 = 4;     %---单元中结点局部编号2对应的结点全局编号列
PosSEleNum3 = 5;     %---单元中结点局部编号3对应的结点全局编号列
PosLEleNum1 = 3;     %---线单元中结点局部编号1对应的结点全局编号列
PosLEleNum2 = 4;     %---线单元中结点局部编号2对应的结点全局编号列

region.cablein = 1;  %---内部电缆实体编号
region.cableout = 3;  %---外部电缆实体编号
region.air = 2;  %---空气实体编号
region.cableinface = 1;
region.cableoutface = [4 7];

% 获取网格基本信息
NumOfcoord = coord(end, PosOfcoord)-coord(1, PosOfcoord)+1;  %---结点数目
NumOfVelement = element_v(end, PosOfelement) - element_s(1, PosOfelement)+1; %---体单元数目
NumOfSelement = element_s(end, PosOfelement) - element_s(1, PosOfelement)+1; %---面单元数目
NumOfLentity = element_l(end, PosOfLentity) - element_l(1, PosOfLentity)+1; %---线实体数目
% 区域单元确定
Numid.air = [find(element_v(:, PosOfVentity)==region.air)];  %---空气域单元编号
Numid.cablein = [find(element_v(:, PosOfVentity)==region.cablein)];
Numid.cableout = [find(element_v(:, PosOfVentity)==region.cableout)];
Numid.cableinface = [find(element_s(:, PosOfSentity)==region.cableinface)];
Numid.cableoutface = [find(element_s(:, PosOfSentity)==region.cableoutface(1));
                      find(element_s(:, PosOfSentity)==region.cableoutface(2))];
% tetramesh(element_v(:, 3:end), coord(:, 2:4), 'FaceColor', [1 1 1]);
% axis equal;

%%
% 单元形状函数以及单元面积
K = zeros(NumOfcoord, NumOfcoord);  %初始化刚度矩阵
b = zeros(NumOfcoord, 1);  %初始化右端项
% volume = zeros(NumOfcoord, 1);

%%
% 确定刚度矩阵K以及右端项b
for ii = 1:NumOfVelement
    volume = 0.5*det(...
        [1 coord(element_v(ii, PosVEleNum1), PosOfcoordx) coord(element_v(ii, PosVEleNum1), PosOfcoordy) coord(element_v(ii, PosVEleNum1), PosOfcoordz); ...
         1 coord(element_v(ii, PosVEleNum2), PosOfcoordx) coord(element_v(ii, PosVEleNum2), PosOfcoordy) coord(element_v(ii, PosVEleNum2), PosOfcoordz); ...
         1 coord(element_v(ii, PosVEleNum3), PosOfcoordx) coord(element_v(ii, PosVEleNum3), PosOfcoordy) coord(element_v(ii, PosVEleNum3), PosOfcoordz); ...
         1 coord(element_v(ii, PosVEleNum4), PosOfcoordx) coord(element_v(ii, PosVEleNum4), PosOfcoordy) coord(element_v(ii, PosVEleNum4), PosOfcoordz)] );  %--计算每个单元的体积

    be = zeros(1, 4);
    ce = zeros(1, 4);
    de = zeros(1, 4);
    sign = 1;
    for jj = 0:3
        be_tmp = [1, 1, 1;
                  coord(element_v(ii, PosVEleNum1+rem(jj+1, 4)), PosOfcoordy), coord(element_v(ii, PosVEleNum1++rem(jj+2, 4)), PosOfcoordy), coord(element_v(ii, PosVEleNum1+rem(jj+3, 4)), PosOfcoordy)
                  coord(element_v(ii, PosVEleNum1+rem(jj+1, 4)), PosOfcoordz), coord(element_v(ii, PosVEleNum1++rem(jj+2, 4)), PosOfcoordz), coord(element_v(ii, PosVEleNum1+rem(jj+3, 4)), PosOfcoordz) ];
        be_tmpValue = det(be_tmp);
        be(jj+1) = sign*be_tmpValue;
        sign = -sign;
    end
    sign = 1;
    for jj = 0:3
        ce_tmp = [1, 1, 1;
                  coord(element_v(ii, PosVEleNum1+rem(jj+1, 4)), PosOfcoordx), coord(element_v(ii, PosVEleNum1++rem(jj+2, 4)), PosOfcoordx), coord(element_v(ii, PosVEleNum1+rem(jj+3, 4)), PosOfcoordx)
                  coord(element_v(ii, PosVEleNum1+rem(jj+1, 4)), PosOfcoordz), coord(element_v(ii, PosVEleNum1++rem(jj+2, 4)), PosOfcoordz), coord(element_v(ii, PosVEleNum1+rem(jj+3, 4)), PosOfcoordz) ];
        ce_tmpValue = det(ce_tmp);
        ce(jj+1) = sign*ce_tmpValue;  
        sign = -sign;
    end
    sign = 1;
    for jj = 0:3
        de_tmp = [1, 1, 1;
                  coord(element_v(ii, PosVEleNum1+rem(jj+1, 4)), PosOfcoordx), coord(element_v(ii, PosVEleNum1++rem(jj+2, 4)), PosOfcoordx), coord(element_v(ii, PosVEleNum1+rem(jj+3, 4)), PosOfcoordx)
                  coord(element_v(ii, PosVEleNum1+rem(jj+1, 4)), PosOfcoordy), coord(element_v(ii, PosVEleNum1++rem(jj+2, 4)), PosOfcoordy), coord(element_v(ii, PosVEleNum1+rem(jj+3, 4)), PosOfcoordy) ];
        de_tmpValue = det(de_tmp);
        de(jj+1) = sign*de_tmpValue;
        sign = -sign;
    end
    

    if(ii<=Numid.air(end) && ii>=Numid.air(1))
        epsilon_this = epsilon.air;
        rho_this = 0;
    elseif(ii<=Numid.cablein(end) && ii>=Numid.cablein(1))
        epsilon_this = epsilon.cablein;
        rho_this = 0;
    elseif(ii<=Numid.cableout(end) && ii>=Numid.cableout(1))
        epsilon_this = epsilon.cableout;
        rho_this = 0;
    end
    K(element_v(ii, PosVEleNum1:PosVEleNum4), element_v(ii, PosVEleNum1:PosVEleNum4)) = K(element_v(ii, PosVEleNum1:PosVEleNum4), element_v(ii, PosVEleNum1:PosVEleNum4))+...
                                                         epsilon_this*(be'*be+ce'*ce+de'*de) / (36*volume);
    b(element_v(ii, PosVEleNum1:PosVEleNum4), 1) = b(element_v(ii, PosVEleNum1:PosVEleNum4), 1)+volume*rho_this / 4;
end
%%
% 添加边界条件
element_cablein = unique( element_v(Numid.cablein, PosVEleNum1:PosVEleNum4) );
K(element_cablein, element_cablein) = 1e20*eye(length(element_cablein));
b(element_cablein) = voltage_cablein*1e20;
element_cableout = unique( element_v(Numid.cableout, PosVEleNum1:PosVEleNum4) );
K(element_cableout, element_cableout) = 1e20*eye(length(element_cableout));
b(element_cableout) = voltage_cableout*1e20;

%%
% 求解
varPhi = K\b;
%%
% 结果云图
% 定义圆形区域的边界
theta = linspace(0, 2*pi, 100);
xc = 1.2*cos(theta)./1e3;
yc = 1.2*sin(theta)./1e3;

% 定义插值网格
[Xq, Yq] = meshgrid(-1.2:0.01:1.2, -1.2:0.01:1.2);
Xq = Xq./1e3;
Yq = Yq./1e3;
Zq = 0.000*ones(size(Xq, 1), size(Xq, 2));

% 使用 inpolygon 确定圆形区域内的点
in = inpolygon(Xq, Yq, xc, yc);

% 使用 griddata 进行插值
Vq = griddata(coord(:, PosOfcoordx), coord(:, PosOfcoordy), coord(:, PosOfcoordz), varPhi, Xq, Yq, Zq);

% 将圆形区域外的值设置为 NaN
Vq(~in) = NaN;

% 绘制云图
pcolor(Xq, Yq, Vq);
axis equal;  % 设置坐标轴比例相等
shading interp;  % 使用插值着色
colormap(jet);  % 设置颜色映射
colorbar;  % 显示颜色条
c = colorbar; % 获取颜色条对象
c.Label.String = 'V'; % 设置颜色条的标签
axis image;
xlabel("m");
ylabel("m");

figure;
linex = linspace(0.5, 1, 20)'./1e3;
liney = zeros(length(linex), 1);
linez = zeros(length(linex), 1);
lineVq = griddata(coord(:, PosOfcoordx), coord(:, PosOfcoordy), coord(:, PosOfcoordz), varPhi, linex, liney, linez);
plot(linex, lineVq, '-*', LineWidth=1);
xlabel("m");
ylabel("V");
% hold on;
% dataExp = readtable("Calculator Expressions Plot 1.csv", "VariableNamingRule","preserve");
% plot((dataExp.("Distance [um]")+500)./1e6, dataExp.("Voltage [V]"));
